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ABSTRACT 

A  general  theory  of  recursive  estimation  is  applied  to  the  particular  problem  of 
estimating  the  time  history  of  a  re-entry  vehicle's  position,  velocity  and  ballistic 
coefficient  from  radar  measurements  of  range,  elevation,  azimuth  and  range  rate. 
The  resulting  algorithms  have  the  following  properties: 

1)  Conceptually  simple  and  easy  to  implement  on  either 
general  purpose  or  special  purpose  computers 

2)  Noniterative 

31  Storage  requirements  independent  of  amount  of  data 

4)  Extremely  fast;  capable  of  real  time  operation 

5)  A  finite  memory  span;  that  is,  can  act  like  a  "sliding 
arc"  algorithm. 

Explicit  formulae  are  provided  along  with  general  discussions  on  the  basic  concept. 


Accept(>d  for  the  Air  P'orec 
Stanley  J.  Wisniewski 
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ALGORITHMS  FOR  ESTIMATING  A  RE-EivlTRY  BODY'S  POSITION, 
VELOCITY  AND  BALLISTIC  COEFFICIENT  IN  REAL  TIME 
OR  FOR  POST  FLIGHT  ANALYSIS 


1.  INTRODUCTION 

This  report  presents  explicit  computational  algorithms  for  reducing  radar 
observations  of  range,  azimuth,  elevation  and  range  rate  to  an  estimate  of  the  time 
behavior  of  a  re-entry  vehicle’s  position,  velocity  and  ballistic  coefficient  (or 
weight-to-drag  ratio,  p).  The  algorithms  are  conceptually  simple  and  computationally 
very  fast.  Computer  storage  requirements  are  small  and  independent  of  the  amount 
of  data  processed.  Iterative  techniques  are  not  used.  The  algorithms  can  employ  a 
finite  effective  memory  span.  Thus  knowledge  of  the  exact  equations  of  motion  and 
an  assumption  of  constant  /3  are  not  required  as  slow  variations  in  )5  can  be  tracked. 
The  algorithms  are  for  real  time  data  processing.  Their  high  speed  also  makes  them 
desirable  for  post -flight  analysis  when  large  amounts  of  data  are  to  be  processed. 
Implemented  on  appropriate  computer  hardware,  they  enable  fruitful  man-machine 
interaction. 

The  algorithms  are  recursive  in  nature  and  the  basic  concept  is  not  new. 

Reference  1  states  it;  Refs.  2  and  3  analyze  its  behavior  for  certain  situations; 

Ref.  4  contains  a  tutorial  discussion  on  the  linear  theory  which  underlies  it;  and  these 
references  are  only  samples;  see  for  example.  Refs.  5  and  6.  The  basic  concept  is 
presently  being  used  for  estimating  position  and  velocity  from  range,  elevation  and 
azimuth  observations  in  the  real  time  system  for  the  Tradex  radar  on  Kwajalein.  * 

The  Kwajalein  algorithm  employs  a  useful  type  of  approximation  which  we  also  discuss. 
The  purpose  of  the  present  report  is  simply  the  collation  of  the  existing  material  into 
a  form  more  directly  suitable  for  implementation.  Hopefully  our  presentation  combines 
formula  and  concepts  almost  ready  for  the  programmer's  pencil  with  generalities 

*  This  application  of  the  basic  concept  was  developed  by  Ken  Ralston  of  Lincoln  Lab. 
independently  from  the  studies  reported  in  Refs.  1-4.  His  work  provided  the 
motivation  for  this  report. 
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sufficient  for  understanding  both  the  concept’ s  power  and  its  shortcomings. 

Notational  requirements  cause  many  of  our  equations  to  appear  rather  foreboding. 
Therefore  we  preview  the  simplicity  of  the  resulting  algorithms  by  giving  Fig.  1. 1, 
which  is  a  verbal  flo  w  diagram  for  one  of  our  algorithms.  (The  same  diagram 
expressed  in  mathematical  symbols  is  given  in  Sec.  5,  Fig.  5. 1.)  The  weighting 
matrices  are  defined  by  recursive  equations  capable  of  real  time  numerical  solution. 
In  a  special  case  of  practical  importance,  explicit,  closed  form  expressions  for  the 
weighting  matrices  are  given. 

An  obvious  and  important  question  is;  how  does  the  accuracy  of  a  fast,  simple 
scheme  such  as  Fig.  1. 1  compare  with  say  an  iterative  solution  of  the  Maximum 
Likelihood  equations?  In  this  report,  we  do  not  address  this  question  as  it  requires 
the  fairly  advanced  mathematics  used  in  Refs.  2  and  3.  However  some  answers  are 
known.  Under  certain  conditions  which  often  occur  in  practice,  our  algorithms  are 
as  efficient  as  Maximum  Likelihood.  These  conditions  are  listed  in  Sec,  3.  Under 
more  general  conditions  (albeit  Gaussian  errors),  the  iterative  Maximum  Likelihood 
approach  intuitively  seems  to  be  better,  but  mathematical  proofs  of  such  intuition  are 
not  available. 

In  Sec.  2  we  go  through  some  preliminaries  and  in  Sec.  3  present  the  two  basic 
algorithms.  Approximations  for  a  special  case  are  considered  in  Sec.  4.  Various 
computational  aspects  are  covered  in  Sec.  5  with  a  final  discussion  given  in  Sec.  6. 


Fig.  1. 1  Flow  diagram  for  one  algorithm 
(see  also  Fig.  5. 1)  . 
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2.  PRELIMNARIES 


We  cryptically  present  some  definitions  and  formulas  which  will  be  used  in  the 
later  developments. 

2. 1.  Basic  Model 


For  the  present,  we  assume  the  re-entry  body  is  moving  in  a  known  force  field  on 
a  trajectory  that  is  completely  specified  by  a  parameter  vector  consisting  of  seven 
elements;  the  six  components  of  position  and  velocity  at  some  instant  of  time  and  P, 
the  weight  to  drag  ratio.  Let  the  vector  p(t) 


P(t)  ■= 


Pl(t) 

PjW 

Po(t) 

P4W 

P=(t) 

p:« 


denote  the  body’s  position  and  velocity  in  a  cartesian,  inertial,  earth  centered 
coordinate  system  where  Pj,(t),  k  =  1,2,3,  are  the  position  components  and  Pj^(t), 
k  =  4, 5,  6,  are  the  velocir  components.  The  dme  behavior  of  this  vector  is  governed 
by  the  vector  system  of  differential  equations 


i  £(t)  =f[E,(t).  t] 

where  the  elements  of  the  vector  ^  are  given  by 


(2.1) 
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iy.  [  *^3  =  Pjc+s^^^  k  =  1, 2, 3 


(2.2) 


fj,Ip(t),tl  = 


,  2,  ,  2,^  1,  3/2 

JPj(t)+P2(t)+P3(t)l 


+ 


d,(t) 


k  =  4,5,6 


where  fi  is  the  gravitational  constant  and  is  the  acceleration  due  to  drag 

given  by 

rsj  1 1  *y 

dj^(t)  Pj^(t)[p^(t)  +  Pg(t)  +  p“(t)]  ^  p(t) 

2W/C^A 


where  /S  =  W/C^ A,  p(t)  is  the  atmospheric  density  and  Pj,(t),  k  =  4, 5, 6  are  the 
inertial  velocities  corrected  to  account  for  the  atmospheric  rotation. 

Define  the  object’s  position  and  velocity  at  time  t  relative  to  a  radar  on  the 
earth’s  surface  by: 


r(t) 

range 

0(t) 

azimuth 

elevation 

r(t) 

range  rate 

m 

azimuth  rate 

m 

elevation  rate. 

Further,  define  the  seven  dimensional  vector  w(t)  by 
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w(t)  = 


^(t) 

Wj(t) 

Bit) 

W,(t) 

m 

w'ct) 

r(t) 

= 

W^(t) 

0(t) 

w_ttl 
O'  ' 

i(t) 

W5(t) 

a(t) 

w,(t) 

where 


ff(t)  = 


P(t) 

P 


Knowledge  of  w(t)  at  t  =  t^  combined  with  Eq.  (2.  i)  enables  the  calculation  of  w(t^^j) 
for  any  value  of  t  , .  Define 


=  (2.3) 

as  the  corresponding  function.  The  evaluation  of  Eq.  (2.3),can  be  done  by: 

•  • 

1)  Conversion  of  r(t  ),d(t  ),  r(t  ),9(t  ),<P(t  )  to  inertial  coordinates 

n  n  n  n  n  n 

eV 

2)  Numerical  integration  of  Eq.  (2. 1)  from  to  t^^^  to  get 

31  Conversion  of  p(t  ,)  back  to  radar  coordinates 

—  n+l 

4)  Evaluation  of  from  a  model  atmosphere. 

The  reasons  for  our  choice  of  the  coordinate  system  w(t)  and  in  particular,  the 

definition  of  Q!(t)  will  become  evident  later  on.  We  discuss  in  Sec.  5  the  use  of  other 

coordinate  systems  and  methods  of  calculating  Eq.  (2.3). 

Define  6  [  t  ,  t  ,w^(t  )]  as  the  7  by  7  matrix  formed  from  the  partial  derivatives 
n  m  “  m 

of  the  elements  of  w(t  )  with  respect  to  the  elements  of  w(t  ),  evaluated  for  w(t  )  = 

—  n  —  m  —  m 
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w  (t  ).  Thus, 
-  m 


Qft  .t  w"(t  )]  = 
n  m  “  m 


(2.4) 


Let  the  q -dimensional  vector  Y(t^)  denote  the  observation  made  by  the  radar  at 

time  t  .  Assume 
n 


Y(t  )  =  H w(t  )  +  v(t  )  (2.5) 

n  “  n  •-  n 

where  H  is  a  qxy  m.atrix  and  v(t  ),  n  =  1, 2, . . . ,  is  a  sequence  of  zero  mean  vector 

—  n 

Gaussian  random  variables  with 


0 


n/^m 


E[v(t  )v'(t  )1  =  < 
—  n  —  m  ' 


Q(t^)  n  =  m 


where  the  prime  denotes  transpose.  To  illustrate  the  use  of  this  model  considei  the 

case  where  at  times  t  ,  n  =  1,  2, . . . ,  the  radar  measures  range,  elevation,  azimuth 

n 

2  2  2  2 

and  range  rate  with  stationary^  independent,  errors  of  variations  cr^,  cr^,  cr^, 
respectively.  Then 


q  =  4 


7 


1 


0  0  0  0  0  0 
0  1  0  0  0  0  0 
0  0  1  0  0  0  0 
0  0  0  1  0  0  0 


Q(g=Q= 


0  0  0 

0  0 

«  -5 ", 

0  0  cr. 


(2.6) 


(2.7) 


Correlation  between  the  errors  (say  between  range  and  range  rate  errors)  are  handled 
by  "fUling  in"  Q(t^). 

Define 


w(t^/m2»  nij^):  The  estimate  of  w(t^)  made  from 

k=mj, . . . ,  m^ 

w(t^/m):  The  estimate  of  w(t^)  made  from 

k  =  l, . ,m  • 

The  purpose  of  this  report  is  the  presentation  of  computational  algorithms  for 

w(t  /m„,m,)  and  w(t  /m).  An  estimate  of  a(t  )  of  course  provides  an  estimate  of 
“  n  2  1  —  n  n 

/3  and  we  often  refer  to  a(t  )  as  an  estimate  of  p. 

n 

2. 2.  Estimation  for  Linear  Dynamical  System 

We  now  change  the  subject  and  consider  some  results  from  the  theory  of  parameter 
estimation  for  linear  dynamical  systems.  Consider  the  system 

x(t  )  =  $(t  ,t  )x(t  )  (2.8) 

—  n  n  n-1  —  n-1  '  ’ 
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where  $  is  a  square,  invertable  matrix.  Assume  we  make  the  observation  2(t^)  of 
the  form 


z(t  )  =  Hx(t  )  +  v(t  ) 

—  n  —  n  —  n 

where  H  and  v(t  )  are  as  defined  and  discussed  with  respect  to  Eq.  (2. 5).  Define 
—  n 

x(t^/m2,  m^)  and  x(t^/m)  in  direct  analogy  with  the  w(t^/m2,mj^)  and  w(t^/m). 

For  the  growing  memory  case  we  have  the  following  formula:* 

x(tyn)=x(t^/n-l)  +  j"^(n/n)H'Q“^(t^) [z(t^)  -Hx(tyn-1)]  (2. 9) 


x(t /n-1)  $(t  .  t  )x(t  ,/n-l) 

—  n  n  n~i  —  n"i 

r(n/n)  =  t^.j)  (2. 10) 

+  H'Q'*(t„)H 

where  "  -1 "  denotes  matrix  inversion. 

For  the  finite  memory  case  we  have  the  following  formula:^ 


x(t  /n.n-r)  =  x(t  /n"l,n-T-l)  + 

-  n  —  n 

i  \n/n,n-T)  {H'Q  ^(t^)[z(t^) -Hx(tyn-l,n-T-l]  - 

$”^\t,t  ,)H'Q‘\t  ,)[z(t  - 

'  n  n-T-i''  ^  '  n-T-r  n-T-r 


*  These  are  explicitly  given  in  Section  3.3  of  Ref.  4. 

t  These  are  not  explicitly  derived  in  Ref.  4  but  the  reference  contains  the  ideas 
necessary  for  their  derivation. 
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where 


r(n/n,n-T)  =  #  ^  (t^,  J(n-l/n-l,n-T-l)  $  \t^)H 

.t  ,)H^Q~V  ,t  ,) 

'  n  n-r-V  ^  ^  n-r-r  '  n  n-r-r 


(2. 12) 


x(tyn-l,n-T-l)  =  #(t^,  t^_^)x(t^_^/n-l,n-T-l) 


The  preceding  formulae  appear  complex  but  ±ey  actually  have  a  simple 
structure  which  we  now  unveil  in  hopes  of  providing  the  reader  with  some  insight.  * 
Consider  first  the  growing  memory  case  of  Eqs.  (2.  9)  and  (2. 10).  Define 


J(n/m): 


I(n): 


The  information  matrix  which 
measures  the  amount  of  infor¬ 
mation  on  x(t  )  contained  in 
z(tj^),  k  =  1, . . . ,  m. 

The  information  matrix  which 
measures  the  amount  of  infor¬ 
mation  on  x(t  )  contained  in 
V  ~  n 


r(n)  =  H’Q'^t^)H 


With  these  definitions,  Eq.  (2. 10)  is  simply  a  statement  that  information  matrices  are 
additive;  i.  e. ,  Eq.  (2. 10)  can  be  rewritten  as 

*  Equations  (2. 9)  through  (2. 12)  are  written  in  a  form  which  leads  most  readily  to  the 
nonlinear  problem  of  real  interest. 
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J(n/n)  =  7(n/n-l)+J(n)  • 


Now  define 


5(n/m):  The  information  vector  which 

is  the  actual  information  on 
x(tn)  contained  in 

k  =  1, ...» m. 

5(n):  Tlie  information  vector  which 

is  the  actual  information  on 

x(t  )  contained  in  z(t  ). 

—  n  —  n 

5(n)  =  H'Q‘\t  )z(t  ) 

—  n  —  n 

B  (n/m)=I(n/m)x(n/m) 

5(n/m)  =  #  ^  (n-l/m) 

x(n/m)=  #(t  ,  t  )x(n-l/m)  . 

—  n  n~l  ~ 


With  these  definitions,  algebraic  manipulation  shows  that  Eq.  (2.  9)  is  simply  a  state¬ 
ment  that  the  information  vectors  themselves  are  also  additive;  that  is,  Eq.  (2.9) 
can  be  rewritten  as 


5(n/n)  =5(iv'n-l)  +5(n)  . 


The  same  type  of  interpretation  can  be  given  to  the  finite  memory  case  of  Eqs, 
(2.  l^and(2. 12).  Define 


/(n/m^,  m^): 


The  information  matrix  for 
z(t, ),  k  =  m,, . . . ,  m- 
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Note  that  J(n/m,  m)  is  the  amount  of  information  on  x(t  )  contained  in  z(t  )  and  thus 

~  n  —  m 


J(n/n,  n)  =  J(n)  . 


Then  Eq.  (2. 12)  can  be  rewritten  as 


7(n/n,  n-T)  =  J(i^n-1,  n-1  -r) 


+  J(n)  -  J(n/n-T-l,  n-T-1) 


that  is,  Eq.  (2. 12)  simply  states  that  we  add  the  information  matrix  for  and 

subtract  the  information  matrix  for  z(t  ,).  Define 

—  n-T-i 

B(n/in^,mj^):  The  information  vector 

for  k  =  mj^, . . . ,  m^ 


Note 


B(n/n,  n)  =  5(n) 


Then  by  algebraic  manipulation,  Eq.  (2. 11)  can  be  rewritten  as: 


B(n/n,n-r)  =5(n/n-l,  n-r-l) 

+5(n)  -5(n/n-T-l,  n-T-1) 


which  says  we  simply  add  and  subtract  the  appropriate  information  vectors. 

The  terms,  information  matrix  and  information  vector,  a^e  used  as  the 
corresponding  concepts  can  be  directly  related  to  information  theory.  *  Reference  7 

*  Our  information  matrix  is  often  called  the  Fisher  Information  Matrix.  The  term, 
information  vector,  is  not  standard  jargon.  It  is  closely  related  to  the  concept  of  a 
sufficient  statistic. 
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IS 


contains  complete  mathematical  discussions  on  this  subject  while  the  Appendix  of 
Ref.  8  contains  a  simplified  discussion. 

The  information  matrix  has  the  following  very  useful  physical  interpretation, 

I  ^(n/m)  =  E  {[x(t  /m)  -x(t  )]  [x(t^/m)  -x(t  )] '  }  (2. 13) 

—  n  — n  —  11  —  n 

that  is,  I  ^(n/m)  is  the  covariance  matrix  of  the  errors  in  x(t^/m).  Reference  4 
discusses  why  the  algorithms  of  Eqs.  (2.9)  and  (2. 10)  give  the  minimum  possible 
I  ^(n/m)  under  the  constraint 


E[x(t  /m))  =x(t  ) . 

—  n  —  n 

An  analogous  interpretation  can  also  be  given  I  ^(n/m^.  ni^^). 
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3.  THE  ALGORITHMS 


We  now  present  the  algorithms  of  interest.  We  do  not  derive  them  in  any 
mathematical  sense,  but  the  motivation  behind  their  choice  is  straightforward. 
Consider  Eq.  (2. 3).  A  Taylor  series  expansion  about  some  estimate  w(t^  gives 


w(t  )=^'[w(t  ),  t  ,t  ] 

—  n  —  n-l  n  n-l 


(3. 1) 


+  Remainder  terms 


where  0  is  the  matrix  of  partial  derivatives  of  Eq.  (2. 4).  Now  assume  w(t^  "  w(t^ 
is  small  enough  to  make  the  remainder  terms  negligible.  It  is  then  reasonable  to  use 
the  linear  data  processing  al^rithms  of  Sec.  2. 2  where  we  make  the  following 
associations: 


x(t  )  ~w(t  )  -^[w(t  ,),t  ,t  ,] 
—  n  —  n  —  n-l  n  n-l 


In  Sec.  2. 1,  we  assumed  that  /3  is  constant  and  that  the  equations  of  motion  of 
the  re-entry  body  are  as  given  by  Eq.  (2.3).  During  re-entry,  this  is  often  a  naive 
assumption.  However,  over  short  periods  of  time,  /5  is  essentially  constant*  and 
the  body's  motion  is  well  approximated  by  Eq.  (2. 3).  Therefore  we  handle  the  real 
problem  of  interest  by  limiting  the  memory  span  of  the  algorithms.  This  gives  a 
tracking  action  whereby  we  can  follow  slow  variations,  in  say 


*  For  the  present  algorithms,  we  do  not  attempt  to  estimate  variations  in  drag  due  to 
body  oscillations.  Thus,  with  respect  to  einy  such  oscillatory  effects,  we  estimate  an 
average  /3. 
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One  basic  algorithm  is  as  follows: 


w(t  /n)  =w(t  /n-l)+J  ^(n/n)H'Q  ^(t  )  [y(t  )-Hw(t  /n-1)] 
—  n  —  n  n-^n~n 


(3.2) 


w(t  /n-1)  =’i'[w(t  /n-l),t  ,t  ] 

—  n  —  n-i  n  n-1 


(3.3) 


J(n/n)=  ,t  ,w(t  /n-l)lX(n-l/n-l)e‘V„>t  ,,w(t  ,/n-l)] 

n  n-1  —  n-1  n  n-i  —  n-i 


+  H'Q  (gH 


for  n  ^  n 


(3.4) 


=  Kn/n) 
o  o 


for  n  >  n 

o 


This  algorithm  is  of  course  analogous  to  the  growing  memory  formulae  of  Eqs.  (2. 9) 

and  (2. 10).  The  main  difference  is  the  "clamping"  of  J(n/n)  at  a  fixed  value , 

J(n  /n  )  for  n  >  n  .  M  J(n/n)  is  not  clamped,  it  goes  to  zero  and  x(t  /n)  converges 
0  0  0  —  n 

to  a  single  unique  trajectory  as  n  increases.  This  would  be  desired  if  ^  were 
constant  and  Eq.  (2. 3)  were  exact.  *  The  clamping  action  limits  the  effective  memory 
span  of  the  algorithm  and  thereby  provides  the  desired  tracking  ability. 

The  second  basic  algorithm  is  directly  analogous  to  Eqs.  (2. 11)  and  (2. 12)  and 
thus  achieves  a  finite  memory  span  without  subterfuge. 


w(t  /n,  n-T)=w(t  /n-l,n-T-l)  + 
—  n  —  n 


J  ^(n/n,  n-T){H' Q  ^(g  [^(g  ”  H w(gn-l,  n-r-l)]  - 


(3.5) 


,)H'Q‘\t  )lx(t  ^  ,)  -Hw(t  /n-l,n-T-l)]} 
n  n-T-l  n-T-i  ^  n-r-l  —  n-r-i 


*  Reference  2  analyzes  the  behavior  of  w(t  /n)  under  these  conditions  for  the  scalar 

”  n 

case  and  gives  conditions  which  guarantee  convergence  and  asymptotic  efficiency. 


(3.6) 


/(n/n, n-r)  =  0'  ^(t  ,t  J J(n-l/n-l,  n-r-i)  0  ^(t  ,t  )  + 

n  n-i  n  n“i 


H'Q'\t  )H  -  e''\t  ,t  )He'\t  ,t  ) 

n  n  n-T-i  n-T-1  n  n-T-i 


w(t  /n-l,n-T-l)  =  ,/n-l,n-T-l), t  ,t  J 

—  n  —  n-l  n  n-i 


(3. 1) 


(3.8) 


where  in  Eqs.  (3. 5)  and  (3.6)  we  have  simplified  the  notation  by  writing 


n  m  n  m  —  m 


Equation  (3. 8)  requires  the  numerical  integration  of  Eq.  (2. 1)  from  t^_j^  backwards  in 
time  to  t^  ^  This  may  be  acceptable  but  if  not,  the  following  formulae  can  be  used 
instead: 

w(t  , /n-l,  n-T-1)  =  w(t  ,/n-2,  n-r-2) 

—  n-T-1  —  n-T-i 

(3.9) 

®  i/n-l, n-T-1)  -w(t  /n-2,n-T-2)]  . 

n-i  n-T-1  —  n-i  —  n-i 


The  algorithm  requires  storage  of  the  observations,  ^(t^)  ;  k  =  n, . . . ,  n-T-1  (and  the 

— (*k  T  ^  used).  However  for  most  cases  of  interest,  this 

requirement  does  not  limit  the  technique's  usefulness. 

As  stated,  both  algorithms  evaluate  ^  each  time  a  new  vector  observation  is 

obtained;  that  is,  we  relinearize  as  per  Eq,  (3. 1)  with  each  new  observation.  However, 

ttl 

in  many  high  data  rate  radars,  relinearization  only  every  r  vector  observation  is 
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required.  In  such  cases  we  can  effectively  combine  the  r,  q-dimensional  observa¬ 
tions  into  a  single  7-dimensional  observation  and  the  algorithm  corresponding  to 
Eqs.  (3. 2)  and  (3. 4)  is 

n 

w(tyn-r)  +  I  \n/n)  ^  {O'  ^(t^,  t^^)  H  'Q  "  H  w(tj,/n-r)j } 

k=n-r+l 

(3.  iO) 


i(n/n)  =0'  ^(t  ,t  )7(n-r/n-r)0  ^(t  ,t  )  + 
n  n-r  '  ' n  n-r 

n  (3. 11) 

k=n-r+l 


where  we  have  again  simplified  the  0  notation.  Equations  (3. 5)  through  (3. 9)  for  the 
finite  memory  case  can  be  modified  in  a  similar  fashion.  This  reduction  in  the  number 
of  relinearizations  can  save  computer  time  as  will  be  discussed  in  Sec.  5. 

Both  algorithms  are  recursive  in  nature.  Initial  conditions  for  the  finite  memory 
algorithm  of  Eqs.  (3.  5)  through  (3. 9)  can  be  obtained  from  the  growing  memory 
formulae  for  n  =  T;  that  is,  we  can  use  Eqs.  (3.2)  through  (3.4)  (for  n^  ^  t)  to  start 
up  the  finite  memory  scheme.  The  initial  conditions  required  for  Eqs.  (3.2)  and  (3.4) 

are  w(t  /O)  and  7(1/0);  that  is,  an  a  priori  guess  on  the  value  w(t  )  and  a  measure  of 

^  -1  ^ 
how  good  the  apriori  guess  is.  [7  (1/0)  can  be  considered  the  covariance  matrix  of 

the  errors  in  the  apriori  guess  of  w(t^)  as  per  the  discussion  of  Eq.  (2. 13).]  If 

7(1/0)  =  0,  then  the  apriori  guess  is  assumed  to  have  'Infinite  errors"  associated 

with  it  and  it  will  not  be  weighed  into  any  future  estimates.  However,  if  7(1/0)  =  0, 

7  \n/n)  will  not  exist  until  the  total  number  of  observations,  nq,  (q  is  dimension  of 

the  observation  vector)  equals  or  exceeds  7.  There  are  various  ways  to  handle  this 

problem  such  as 
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1)  Make  J(l/0)  small  but  not  zero 

2)  Combine  the  first  few  observations  into  a  single  vector 
observations  as  per  the  discussions  of  the  preceding 
paragraph. 

The  choice  depends  on  the  particular  problem  at  hand. 

In  Sec.  2.2  we  interpreted  J  ^(n/n)  and  I  ^(n/n.n-r)  as  the  covariance  matrices 

of  the  errors  in  the  estimates,  x(t  /n)  and  x(t  /n,  n-r).  The  errors  in  the  estimates 

^  ^  ^  “1 

w(t^/n,  n-r)  can  be  similarly  interpreted  in  terms  of  I  (n/n,  n-r)  of  Eq.  (3.  6) 

provided: 

1)  The  estimate  errors  are  small  enough  so  that  the  remainder  terms 
of  Eq.  (3. 1)  are  truly  negligible  and 

2)  is  "nearly"  constant  and  Eq.  (2. 3)  is  a  "good"  approximation 
over  the  memory  span. 

The  relationship  between  the  errors  in  w(t^/n)  for  n  ^  n^  and  J  ^(n^/n^)  have  not 
yet  been  investigated  but  under  the  above  two  conditions,  they  should  be  "in  the  same 
general  ball  park. " 

The  two  conditions  of  the  preceding  paragraph  are,  of  course,  the  conditions 
which  lead  to  the  choice  of  the  algorithms  themselves.  When  the  conditions  are  not 
satisfied,  it  is  very  difficult  to  perform  explicit  error  analyses  but  the  algorithm 
performance  should  degrade  gracefully  as  the  conditions  are  violated.  If  the  algorithm 
fails,  it  will  probably  be  due  to  one  or  a  combination  of  the  following  factors: 

1)  An  incomplete  error  model;  i.  e. ,  the  neglect  of  error  sources 
such  as  bias  errors  or  correlations 

2)  Insufficient  data;  i.  e. ,  good  estimates  are  impossible  using  just 
the  observations  made  over  the  memory  span 

3)  Bad  initial  conditions. 

The  recursive  seherne  can  be  modified  to  handle  a  more  extensive  error  model.  The 
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problem  of  insufficient  data  can  be  solved  by  buying  a  more  expensive  radar.  If  the 
choice  of  initial  conditions,  w(tj/0),  is  poor,  the  algorithm  may  search  around  for  a 
while  until  it  locks  on  to  the  trajectory.  If  w(t^/0)  is  very  bad,  it  may  fail  to  lock  on 
at  all.  For  post  flight  analyses  this  is  no  problem  as  different  initial  conditions  can  be 
tried.  *  For  real  time  applications,  failure  to  lock  on  must  be  avoided  and,  if 
necessary,  an  auxialiary  algorithm  used  to  provide  satisfactory  initial  conditions. 

Our  use  of  terms  such  as  "transient, "  "search,  "  "lock  on"  and  "track"  is 
deliberate  as  the  algorithms  are  directly  analogous  to  a  feedback  control  system 
(or  servomechanism)  designed  to  track  an  input  signal  and  its  time  derivatives.  This 
analogy  is  mentioned  as  the  author  has  found  it  to  be  a  valuable  aid. 


*  For  example,  an  initial  transient  can  be  removed  by  feedback  of  the  w(t  /n) 
obtained  after  lock  on. 
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4.  A  SPECIAL  CASE 


In  Sec.  3  we  gave  the  algorithms  in  a  general  form.  We  now  restrict  interest  to 
a  special  case  and  introduce  some  approximations  of  the  type  found  in  the  Kwajalein 
algoridim.  This  enables  certain  closed  form  expressions  which  greatly  reduce 
computation  time  and  complexity.  In  addition,  some  of  the  behavioral  prqjerties  of 
the  algorithms  can  be  displayed. 

Assume  we  observe  range,  elevation,  azimuth  and  range  rate  and  that  the  errors 
in  these  observations  are  independent  of  each  other  and  are  stationary.  Then  we  have 
H  and  Q(tj^)  =  Q  as  given  by  Eq.  (2. 6)  and  Eq.  (2. 7).  The  case  where  the  range  rate 
observation  is  not  present  is  included  as  a  special  case.  Assume  further  that  the 
observations  are  equally  spaced  in  time  by  A  units  of  time.  Thus 


t 

n 


=  t  ,  +  Ak 
n-k 


The  crucial  approximation  is  to  assume  0[t  ,  t  , ,  w(t  )]  is  of  the  form 

n  n  ‘*xC  II  “ic 
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where 


y  =v(t  ,)r(t  J 
r  n-k  n-k 


n-k" 


n-k' 


n-k" 


n-k" 


n-k" 


It  is  easily  seen  that  the  y^,  and  y^  are  simply  the  coefficients  that  split  the 
total  acceleration  due  to  drag  into  its  respective,  r,  9  and  0  coordinates. 

The  motivation  behind  the  choice  of  Eq.  (4. 1)  is  as  follov/s.  Consider  a  small 
variation,  5w(t+Ak),  in  w(t-}-Ak)  due  to  small  variations,  5w(t),  in  w(t).  Then  to  a 
first  approximation 


5  w(t+Ak)  =  0(t+Ak,  t)  6  w(t) 


where  we  have  again  simplified  the  0  notation.  Now  consider  the  range  coordinate. 
Equation  (4. 1)  states. 


5  r(t+Ak)  =  6  r(t)  +  Ak  6  r(t)  + 


(Ak)^ 

2 


y  6a(t) 
r 


that  is  we  have  made  the  approximations 


3r(t+Ak)  ^ 
ar(t) 

9r(t+Ak)  _ 

9t(t)  ■ 

9r(t+Ak)  _  (Ak)^ 

9a  (t)  ■  2  ”^r 


(4.2) 
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and  in  addition,  have  neglected  the  effect  of  the  perturbations,  6  0(t),  5  0(t),  6  ^(t) 
and  6  0(t).  These  perturbation  effects  are  dropped  because  the  corresponding  partial 
derivatives  are  small;  a  valid  step  provided  the  corresponding  perturbations,  5  0(t), . . 
are  not  too  large.  The  incorporation  of  (Ak)  behavior  only  with  respect  to  5Qf(t), 
is  motivated  by  the  assumption  that  in  practice,  5a(t)  will  be  relatively  large. 
Equation  (4. 1)  models  variations  in  constant.  Section  5  discusses  the  use  of 

exponential  atmospheres  and  Eq.  (5. 2)  can  be  used  to  check  this  approximation. 

The  y^,  y^  and  y^  are  time  varying  quantities  but  in  many  applications,  they 
will  be  essentially  constant  over  the  memory  span  of  the  algorithms.  Thus  we  make 
the  additional  assumption  that  they  are  constants,  denoted  by  Tj..y^f  and  let 
0(tn,  t^  denote  the  corresponding  version  of  Eq.  (4. 1).  With  this  final  approxima¬ 
tion,  it  is  possible  to  obtain  closed  form  analytic  expressions  for  J  ^(n/n)  and 
I  ^(n/n,  n-r). 

Under  our  various  assumptions, 

J(n/n,n-T+l)  =  J(T/T)  . 

Thus  we  need  consider  only  X  ^(n/n).  Equation  (3. 4)  can  be  written  as  (for  n:^n^) 

n-1 

k=0 

Evaluation  of  J  \n/n)  proceeds  easiest  by  use  of  the  orthonormal  polynomials  defined 
by 

V  n 
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ial 


of 


n(n  -1) 


] 


‘2  I  kA  - 


j_  /  180  r 

2  !  ">  2 
A  V  n(n*"-l)(n  -4) 


180  I  ,,  ^,2  ,  2  ^2 

(kA)  -(n-l)kA  +A 


where 


e 


a- 


/ 


^  A.(l:A)X.(kA)  =  / 
kM3 


1 

0 


i=j 


These  orthonormal  polynomials  are  discussed  in  many  places,  see, 
Refs.  9  and  10  (which  are  orthonormal  over  k=  1, . . . ,  n).  Define 


) 


ed 


(n-l)(n-2) 

6 


(4.4) 


for  example, 
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D  = 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


(n-l)A 

2 


0 


0 


0 


0 


0 


0 

(n-l)A 

2 

0 

0 


0 


0 


y^A>-l)(n-2) 

12 

rgA^(n-l)(n-2) 

12 


(n-l)A  (”-^Xn-2) 

2  12 


0 


0 


y^(n-l)A 


yg(n-l)A 


y^(n-l)A 


Then 


Eq.  (4.6) 


e"\t  ,t  JD  =  iXkA)  . 
n  n-k 


(4.7) 


Substitution  of  Eq.  (4.  7)  into  (4. 3)  gives 


I  ^(n/n)  =  DI^D' 


(4.8) 


where 


n-l 


1^=^  A'(kA)H'Q'^HA(kA) 


k=0 
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Equations  (2. 6)  and  (2. 7)  for  H  and  Q,  and  the  orthonormal  properties  of  the  A.(kA) 
Eq.  (4.4),  result  in  being  a  diagonal  matrix  with  main  diagonal  elements,  J„, 
given  by 


^11  = 


n 


-22 


n 


^33 


n 


^44  = 


^55  = 


'66 


*77 


2  2 
A  n(n  -1) 

120^ 

r 

2  2 
A  n(n  -I) 

2  2  . 
A  n(n  -I) 

0 

12cr 


+ 


n 

2 

cr. 

r 


<P 

2  4  2  2 

I  A  n(n  -l)(n  -4) 

720 


y^A^n(n^-l) 

O'?  12 

r 


(4.9) 


where 


I 


2 


(4. 10) 


Since  the  inversion  of  diagonal  matrices  is  considered  to  be  trivial,  the  calculation  of 
the  final  closed  form  expression  for  I  ^(n/n)  now  proceeds  directly  by  doing  the 
matrix  multiplication  of  Eq.  (4.8)  using  Eqs.  (4.6)  and  (4.9).  In  order  to  display  the 
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<A), 


9) 


0) 


)f 


result  in  a  relatively  simple  form,  we  partition  J  (n/n)  as 


where 


T  ^(n/n)  = 


A  A  A  g 
PP  pv  vP 

A  A  _ 
pv  w  vp 

A'  A' o  Ago 

P/3  /3/3 


(4.11) 


PP 


vv 


^/3 


and  A  ,A  o>A 
pv’  p^’  vp 


3  X3  matrix  corresponding  to  the  position  coordinates 
3  x3  matrix  corresponding  to  the  velocity  coordinates 
ixi  matrix  corresponding  to  the  Q:(t)  coordinate 

are  the  resulting  off-diagonal  matrices.  Define 


r= 


Then 


2  2 
(n-l)V 


47 


44 


0 


0 


0 


n(n+l) 


0 


0 


+ 


A^(n-l)^(n-2)^ 


1447, 


77 


IT' 


n(n+ 1) 
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Note  that  Eqs.  (2. 6)  and  (2. 7)  give 


Thus  the  4X7  matrix  I  \n/n)H'Q  ^  can  be  calculated  by  inspection  from  Eq.  (4. 11). 
This  4x7  matrix  can  then  be  substituted  directly  in  Eq.  (3.  2)  to  completely  specify 
that  algorithm.  The  necessary  expressions  for  the  algorithm  of  Eq.  (3. 5)  are  obtained 
in  a  similarly  easy  fashion. 

The  case  where  only  range,  azimuth  and  elevation  are  observed  can  be  obtained 

-1  2 

from  our  expressions  for  T  (n/n)  simply  by  letting  o.  —“o.  When  this  is  done, 

-1  ^ 

T  (n/n)  develops  many  more  symmetries. 

The  approximation  of  Eq.  (4. 1)  and  the  constancy  of  the  and  are  valid 

assumptions  over  only  short  periods  of  time.  Thus  if  their  consequences  are  to  be 

used,  the  effective  memory  span  (t  or  n^)  of  the  algorithm  cannot  be  allowed  to 

become  too  large  even  if  /3  is  constant  and  the  equations  of  motion  are  known  over 

longer  periods  of  time. 
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The  approximation  of  Eq.  (4. 1)  is  also  useful  in  the  insight  it  provides  on  the 
performance  of  the  algorithms.  Assume,  as  per  the  discussions  of  Sec.  3,  that 
conditions  are  such  that  J  ^(h/n)  represents  the  covanance  matrix  of  the  errors. 
An  instructive  case  occurs  when  doppler  is  not  present  and  the  re-entry  body's 


velocity  vector  lies  entirely  in  the  radial  direction  to  give  =  0- 


simplicity  we  also  assume  n  »  1.  Then  Eq.  (4. 11)  reduces  to 


For 
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Eq.  (4.12) 


le 


36a  60a 


An  important  question  is  the  performance  of  our  algorithms  if  we  use  them  out¬ 
side  the  atmosphere  where  we  actually  should  not  try  to  estimate  P .  Such  a  procedure 
would  alleviate  the  need  for  an  algorithm  change  as  the  re-entry  occurs.  The  resulting 
degradation  can  be  measured  by  repeating  our  analyses  for  the  case  where  P  need  not 
be  estimated.  The  result  is  as  follows;  if  we  use  our  algorithms  exo-atmospherically, 
we  unnecessarily  degrade  the  range  accuracy  by  ^  9/4  ~  1.5  and  the  range -rate 
accuracy  by  V 192/12  =  4.  The  amount  of  degradation  can  be  somewhat  misleading 
as  we  are  estimating  the  range  and  range  rate  at  the  time  of  the  last  observation  and 
at  this  point,  the  range -range  rate  estimate  errors  are  highly  correlated.  The  picmre 
changes  if  we  compare  the  errors  in  the  estimates  of  range  and  range  rate  as  calculated 
at  the  midpoint  of  the  memory  span;  that  is,  if  for  the  finite  memorj'  algorithm,  we 

consider  the  errors  in  r(t  /n,  n-r)  and  r(t  /n,  n-r)  where  t  =  t  +  At/2.  This 

m  m  m  n-T 

is  the  case  of  interest  in  post  flight  analysis.  At  midpoint,  the  range -range  rate 

estimate  errors  are  uncorrelated;  the  degradation  in  range  is  still  V9/4  »  1. 5; 

but  there  is  no  degradation  in  range  rate  accuracy.  These  results  are  taken  from  Ref. 

9  which  contains  curves  showing  the  effect  of  the  degree  of  polynomial  and  point  of 

estimation  on  the  errors  in  parameter  estimates.  These  curves  are  applicable  in  the 

special  cases  under  discussion.  Note  that  our  comparisons  are  for  the  special  case 

Ta  =  7 ,  =  0.  It  is  seen  from  Eq.  (4. 11),  that  if  instead,  y  =  0,  y  ,  /  0,  y-  /  0,  the 
0  <p  I  (p  y 

estimates  of  range  and  range  rate  are  not  effected  at  all  by  the  p  estimate.  In  this 
sense,  our  comparisons  are  for  the  worst  case. 
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5.  COMPUTATIONAL  ASPECTS 


The  two  basic  algorithms  were  giver  in  Sec.  3.  For  completeness,  we  give  the 
corresponding  flow  diagrams  in  Figs.  5. 1  and  5. 2.  It  is  seen  that  Fig.  5. 1  is  merely 
Fig.  1.1  expressed  in  terms  of  mathematical  symbols  instead  of  verbal  descriptions. 
The  coordinate  conversions  indicated  in  Fig.  1. 1  are  not  explicitly  shown  as  the  w(t) 
coordinate  system  is  based  on  r^dar  coordinates  and  multiplication  by  the  matrix  K 
performs  the  necessary  coordinate  conversions. 
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Old  estimate 


New  observation 


Fig.  5. 1  Flow  diagram  for  Eqs.  (3.2)  -  (3.4)  . 
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Old  estimate  New  observation 

w{t  /n-l,n-T-i)  y(t  ) 

—  n-1  ^  n 


Fig.  5.2  Flow  diagi-am  for  Eqs.  (3.5)  -  (3.8)  . 
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These  algorithms  are  not  quite  ready  for  the  programmer’s  pencil  in  that  certain 
questions  of  implementation  can  only  be  resolved  using  the  particulars  of  the  actual 
problem.  We  now  give  a  conglomerate  of  discussions  relating  to  such  mundane,  but 
very  important  questions. 

We  have  presented  two  different  algorithms.  The  algorithm  of  Eqs.  (3. 2)  through 
(3.4)  is  the  easiest  to  iix'.plement  and  requires  the  least  amount  of  computation  time, 
but  it  is  basically  an  ad  hoc  procedure.  The  algorithm  of  Eqs.  (3. 5)  through  (3. 9)  uses 
a  chosen  finite  memory  span  but  this  advantage  is  obtained  at  the  expense  of  a  more 
complicated  algorithm  which  requires  some  storage  and  an  increased  computation  time 
which  may  or  may  not  be  important.  Cftjviously  the  choice  between  the  two  algorithms 
depends  on  the  problem  of  interest.  However,  for  "scientific"  post  flight  analyses,  the 
preciseness  of  Eqs.  (3. 5)  -  (3. 9)  appears  to  override  any  computational  disadvantages. 

The  parameters  n^  and  t  determine  the  effective  memory  span  of  the  algorithm. 

Let  Tj^  be  the  maximum  time  interval  over  which  ^  is  essentially  constant  and 
Eq.  (2. 3)  is  a  good  approximation.  Let  be  the  maximum  time  interval  over  which 
the  approximations  of  Sec.  4  are  valid.  Let  T  =  T^^^,  or  T  =  min  {Tj^,T2}  depending  on 
whether  or  not  the  equations  of  Sec.  4  are  used.  It  might  be  considered  reasonable  to 
choose  memory  spans  corresponding  to  T.  However,  the  problem  can  rarely  be  so 
easily  resolved.  First  of  all,  T  is  difficult  to  specify  as  models  such  as  a  constant  /3 
or  approximations  a  la  Sec.  4  almost  always  fail  gradually;  that  is,  a  small  increase 
in  T  means  the  model  and  approximations  are  only  slightly  worse.  Secondly,  even  if 
T  were  specified,  a  longer  memory  span  might  still  be  desirable  to  increase  the  ability 
of  the  algorithm  to  smooth  out  the  effects  of  measurement  noise.  Thus  compromises 
between  the  effects  of  model  and  approximation  errors  and  the  effect  of  the  measurement 
errors  might  be  required.  A  third  complication  is  the  presently  unknown  relationship 
between  n^  and  the  effective  memory  span  of  the  algorithm  of  Eqs.  (3. 2)  -  (3. 4).  A 
final  difficulty  is  the  altitude  and  geometry  dependence  of  T.  Thus  even  for  a  given 
re-entry,  no  single  value  of  n^  or  t  need  be  optimum  for  the  entire  flight.  The 
overall  outlook,  however,  is  not  really  so  dismal  as  the  memory  span  is  an  easily 
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varied  parameter  and  good  values  can  be  learned  through  experience.  A  pre- 
programmed  memory  span  dependence  on  altitude  and  slant  range  is  easy  to  provide. 
Furthermore,  in  many  instances  the  ultimate  in  accuracy  is  not  needed  and  overly 
conservative  (i.  e.  small)  memory  spans  are  acceptable.  If  the  ultimate  is  desired 
for  post  flight  analysis,  a  man  can  monitor  the  results  of  repeated  processing  to  see 
which  memory  span  works  best  and  make  the  appropriate  adjustments. 

The  matrix,  Q(t  ),  is  fixed  by  the  observation  error  distributions.  Time 
n 

variations  in  Q  may  arise  in  various  ways;  for  example,  when  the  signal -to -noise 
ratio  is  used  to  measure  the  variances  of  the  observation  errors.  Note,  however, 
that  it  is  the  relative  magnitude  of  the  various  elements  of  Q  that  is  important,  not 
their  absolute  value  (unless  I  ^(n/n)  is  to  be  associated  with  the  estimate  errors). 

For  many  applications,  a  constant  Q  will  be  satisfactory. 

The  algorithms  are  based  on  the  assumption  of  white  observation  errors;  that  is, 
statistical  independence  in  time.  This  is  often  not  true  as,  for  example,  range  and 
angle  tracking  loops  may  be  incorporated  within  the  radar  itself.  If  the  memory  span 
of  the  algorithm  is  sufficiently  long  compared  to  the  error  correlation  time,  little 
harm  is  done  although  I  ^(n/n)  is  then  no  longer  the  covariance  matrix  of  the  errors.  * 
However,  it  is  often  advantageous  to  simply  feed  the  algorithm  with  observations 
taken  far  enough  apart  in  time  to  be  "effectively"  uncorrelated.  If  these  time  intervals 
are  different  for  different  observation  coordinates  such  as  range  and  angle,  the  ideas 
behind  Eq.  (3. 10)  can  be  used.  Alternately,  the  values  of  Q  can  be  modified  to 
correspond  to  an  "equivalent"  white  noise  process. 

Our  algorithms  estimate  the  body's  state,  w(tX  at  time  t^  given  the  observations 
up  to  time  t^.  For  post  flight  analysis,estimates  at  a  time,  t^,  in  the  middle  of  the 
memory  span  are  of  more  interest.  w(t^/n)  can  be  calculated  from  w(t^/n)  by 
numerical  integration  backwards  from  t^  to  t^.  However,  the  algorithms  can  be 
modified  to  provide  midpoint  estimates  directly. 

*  See,  for  example.  Ref.  11  which  discusses  the  asymptotic  efficiency  of  least 
squares  processing  (i.e.  a  white  noise  filter)  in  the  presence  of  correlated  errors. 
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We  have  presented  algorithms  that  estimate  the  re-entry  body’s  position  and 

velocity  in  terms  of  radar  coordinates  and  aft  ),  i.  e.  the  w(t  )  coordinate  system. 

n  —  n 

The  choice  of  radar  coordinates  was  made  for  several  reasons.  No  coordinate 

conversions  of  the  raw  data  are  required.  This  is  important  in  high  data  rate  radars 

when  the  number  of  relinearizations  is  reduced  using  the  ideas  of  Eq.  (3. 10).  It  also 

simplifies  the  incorporation  of  range  rate  and  gives  a  simple  form  for  Q.  Finally, 

the  residuals  between  the  observations  and  the  smoothed  trajectory  are  automatically 

available  and  they  are  important  in  post  flight  analysis.  However,  these  advantages 

might  be  overweighted  when  the  approximations  of  Sec.  4  are  used  as  Eq.  (4. 1)  might 

be  a  better  approximation  in,  say,  an  inertial  coordinate  system  than  in  radar 

coordinates.  If  estimation  in  a  nonradar  coordinate  system  is  desired,  the  algorithms 

can  be  appropriately  modified  by  simply  considering  the  coordinate  converted  radar 

observations  as  the  actual  observations.  The  evaluations  of  the  corresponding  Q  \t  ) 

n 

can  be  done  using  the  well-known  partial  derivatives  of  the  new  coordinate  system 
with  respect  to  the  radar  coordinates.  *  If  range  rate  observations  are  to  be  used,  one 
approach  is  to  do  the  coordinate  conversions  using  an  elevation  and  azimuth  rate 
calculated  from  the  past  estimates  and  then  remove  their  effect  by  assigning  these 
"fake"  observations  infinite  error  variances. 

P(tn) 

We  estimate  ^ —  rather  than  /?  where  p(t^)  is  the  air  density.  This  is 

valuable  at  high  altitudes  where  p(t)  —  0  as  it  keeps  various  "gains"  from  "blowing 
up.  "  It  has  another  advantage  in  that  a  simple  exponential  atmosphere  may  prove 
satisfactory  for  the  actual  processing.  For  example,  if  we  use  the  model 


*  The  Kwajalein  algorithm  estimates  position  and  velocity  in  an  x,  y,  z  system  without 
appropriate  modification  of  the  Q  matrix.  However,  since  range  rate  is  not  used  and 
/3  is  not  estimated,  it  can  be  proved  that  this  does  not  degrade  the  performance  of  the 
algorithm.  Unfortunately  a  proof,  although  straightforward,  is  not  available  in  the 
unclassified  literature. 
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where  p  and  h  are  scale  factors  and  h(t  )  is  the  height  of  the  body  above  the  Earth 
0  0  n 

at  time  t  ,  then 
n 


|-h(g-h(t^-i)-| 

'l  h  J 


Qr(t)  =  a!(t  )e 
n  n-l 


(5.2) 


Because  of  the  tracking  nature  of  the  algorithms,  Eq.  (5. 2)  need  be  valid  only  over 
n^  the  memory  span  of  the  algorithm  and  the  value  of  need  not  be  known.  A  pro¬ 

grammed  variation  of  h^  with  altitude  could  also  be  included.  Of  course  to  estimate 
/3  from  ^  standard  atmosphere  should  be  used.  ('Fhis  use  of  an  exponential 

model  is  especially  fruitful  in  some  real  time  applications  v/here  )3  itself  is  just  a 

nuisance  parameter.)  However,  the  decision  to  estimate  Q!(t  )  is  not  sacred  and  can  be 

n 

modified  to  estimate  /3  or  some  other  function  of  )3  if  the  need  arises.  For  example, 
a  mach  number  dependence  might  be  incorporated. 

We  discussed  the  evaluation  of  Eq.  (2.3)  in  Sec.  2  using  numerical  integration 
in  an  inertial  coordinate  system.  However,  in  practice  the  choi"'  of  this  coordinate 
system  depends  primarily  on  the  programmer's  whims.  Furthermore,  because  of 
the  algorithm's  tracking  action,  the  model  for  Eq.  (2.3)  need  be  valid  only  over  the 
algorithm's  memory  span.  Thus  analytic  approximations  to  these  equations  of  motion 
may  be  acceptable  in  place  of  any  numerical  integration.  This  argument  is  analogous 
i  to  that  used  for  the  exponential  atmosphere. 

The  discussions  of  Sec.  4  are  not  the  last  definitive  word  on  error  analysis  but 
they  indicate  that,  in  many  applications,  it  will  be  satisfactory  to  use  our  algorithms 
outside  the  atmosphere.  However,  if  desired,  we  can  employ  a  preprogrammed 
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algorithm  change  to  begin  estimating  p  at  some  predetermined  altitude.  The  position 
and  velocity  estimates  obtained  exo -atmospherically  furnish  the  initial  conditions  for 
this  re-entry  phase.  The  algorithms  for  estimating  just  position  and  velocity  exo- 
atmospherically  are  an  obvious  special  case. 

In  Sec.  4,  we  presented  an  approximation  to  Q[z  ,t  ,w(t  )] .  Although  this 

n  m  —  m 

approximation  appears  to  have  a  wide  range  of  application,  situations  may  definitely 
occur  where  a  more  accurate  evaluation  is  required.  Tliere  are  several  ways  to 
proceed.  More  accurate  analytic  approximations  can  be  developed.  Exact  analytic 
albeit  complex,  formulae  are  probably  possible  assuming  an  exponential  atmosphere. 
The  most  general  approach  however,  is  the  use  of  numerical  techniques  and  of  the 
various  possibilities,  the  following  appears  to  be  the  most  satisfactory.  In  concept, 
at  least,  we  can  write  the  7-dimsnsional  vector  system  of  first-order,  nonlinear 
differential  equations 


dt =£fw(t),  t] 


where  the  elements  of  g,  g^[  w(t),  t] ,  k  =  1, . . . ,  7  include  dependence  on  the  atmospheric 
density,  earth  rotation,  etc.  Define  the  seven  by  seven  matrix 


9gj[  w(t),  t] 
9^^ 


B[w°(t),t]  = 


9g7[  w(t),  t] 

9^^ 


9gjw(t),t] 

9w^(ty 


9g^[w(t),t] 


where  the  partial  derivatives  are  evaluated  for  w°(t)  as  calculated  from 
we  have  the  matrix  differential  equation 


Then 
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with  the  initial  conditions 


m  m  —  m 


which  can  be  solved  by  numerical  integration  from  t  tot.*0^[t,t  ,  w(t  )]  can 

m  n  n  m  m 

be  found  by  numerical  matrix  inversion  or  by  integrating  backward  from  t^  to  t^. 

For  the  general  cases  not  covered  by  Sec.  4,  the  7(n/n)  matrix  can  be  numerically 
calculated  by  the  recursive  formula,  Eqs.  (3. 4)  or  (3. 6).  Numerical  matrix  inversion 
is  required  to  obtain  I  \n/n).  There  are  various  ways  to  reduce  the  computation 
time  required,  see  for  example  the  discussions  in  Appendix  A.  7  of  Ref.  4.  The 
number  of  matrix  inversions  can  be  reduced  by  relinearizing  only  every  r^^  vector 
observation  as  illustrated  by  Eq.  (3. 10). 

J(n/n)  of  Eq.  (3. 4)  becomes  constant  for  n  >  n^.  Similarly,  in  many  cases, 

J(n/n,  n-r)  becomes  essentially  constant.  In  some  applications,  it  might  be  acceptable 
to  simply  precalculate  these  constant  values  and  just  use  them  in  the  algorithms.  This 
i’ltroduces  additional  transients  into  the  system  but  they  may  not  be  important. 
Similarly,  in  the  special  case  of  Sec,  4  it  might  be  satisfactory  to  merely  use  the  for¬ 
mulae  for  the  case  where  n  »  1  which  simplifies  the  expressions. 


*  In  practice,  the  integration  would  probably  be  done  in  inertial  coordinates  and  then 
converted  into  the  w(t)  system. 
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6.  DISCUSSION 


The  algorithms  of  this  report  have  not  yet  been  programmed  and  used.  However, 
the  basic  underlying  ideas  have  been  proven  to  be  sound  by  analysis,  simulation  and 
field  application.  The  approximation  of  Sec.  4  is  physically  reasonable  and  leads  to 
an  intuitively  satisfying  system.  We  therefore  feel  the  algorithms  are  ready  for 
implementation  using  the  approximations  of  Sec.  4.  Questions  such  as  memory 
length  and  actual  performance  are  best  resolved  by  running  the  algorithms. 

We  have  purposely  addressed  ourselves  to  a  ver>'  special  problem,  namely  the 
estimation  of  seven  parameters,  position,  velocity  and  ^  from  a  single  radar. 

However,  it  is  obvious  that  the  theory  is  actually  applicable  in  a  far  more  general 
context.  By  introducing  some  more  coordinate  conversions  and  allowing  variations  m 
the  matrix  H,  the  technique  can  be  used  equally  well  for  multiple  sites,  each  with 
different  tracking  abilities,*  for  example,  the  combination  of  Baker -Nunn  and  radar  data. 
In  addition,  many  other  parameters  can  be  included  in  the  model  such  as  radar  bias 
errors,  station  location  errors  and  more  complex  parameterizations  of  the  body's 
motion.  These  various  extensions  can  be  implemented  directly  using  the  numerical 
techniques  of  Sec.  5  to  evaluate  ©  and  r,  but  by  some  diligence,  it  is  probable  that 
simplified  analytic  approximations  such  as  those  given  in  Sec.  4  can  also  be  developed. 
In  fact,  extension  to  a  parameterization  of  a  time-varying  $  such  as 

P 

j«l 


where  the 


are  considered  unknown  parameters  requires  a  fairly  trivial  analysis 


provided  the  orthonormal  polynomials  of  Sec.  4  are  employed.  The  same  is  true  for 


the  incorporation  of  coherent  range  acceleration  and  the  estimation  of  lift  forces.  In 


addition,  if  good  dcppler  is  available,  it  presently  appears  reasonable  to  estimate 
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body  oscillations  by  modeling  ^  in  terms  of  a  sinusoid  of  unknown  frequency.  *  The 

tJl 

basic  recursive  schemes  can  also  be  generalized  to  handle  n  order  Markovian  as 
well  as  white  observation  noise. 

A  critical  parameter  of  the  data  processing  algorithms  is  the  memory  span  as 
determined  by  n^  or  r.  We  have  discussed  the  use  of  a  preprogrammed  altitude 
dependence  for  varying  these  quantities  and  also  the  use  of  a  m*- n  to  determine  the 
best  values  in  post  flight  analysis  by  repeated  processing.  However,  the  most 
interesting  and  powerful  approach  would  be  to  use  the  radar  data  itself  to  determine 
the  memory  span.  The  basic  tenets  of  such  a  technique  are  already  known  and  are 
applicable  not  only  to  rnemor^^  span  control  but  also  to  problems  such  as  an  automatic 
transition  from  an  exo -atmospheric  to  a  re-entry  algorithm.  The  analysis  of  the 
technique  will  be  helped  by  the  analogy,  mentioned  in  Sec.  3,  between  the  basic 
recursive  algorithms  and  feedback  control  theory.  As  with  most  nonlinear  control 
systems,  simulation  on  a  computer  complex  capable  of  man -machine  interaction 
will  probably  be  required.  This  may  be  an  area  of  future  research. 


*  Recursive  techniques  similar  to  those  of  this  report  exist  for  estimating  an  unknown 
frequency.  A  Lincoln  Laboratory  report,  ’’Computationally  Feasible  Frequency 
Estimates  in  the  Presence  of  Unknown  Phase  and  Amplitude’  by  L. A.  Gardner,  Jr. , 
is  presently  in  preparation. 
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